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The stability of doubly quantized vortices in dilute Bose-Einstein condensates of 23 Na is exam- 
ined at zero temperature. The eigenmode spectrum of the Bogoliubov equations for a harmonically 
trapped cigar-shaped condensate is computed and it is found that the doubly quantized vortex 
is spectrally unstable towards dissection into two singly quantized vortices. By numerically solv- 
ing the full three-dimensional time-dependent Gross-Pitaevskii equation, it is found that the two 
singly quantized vortices intertwine before decaying. This work provides an interpretation of recent 
experiments [A. E. Leanhardt et al. Phys. Rev. Lett. 89, 190403 (2002)]. 
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The stability of single-quantum vortices has been stud- 
ied extensively after such vortices were first observed in 
dilute alkali atom Bose-Einstein condensates (BECs) . 
Since a single-valued complex order parameter describes 
the state of the condensate, its phase must undergo a 2im 
change along a loop encircling a vortex, where n is the 
quantum number of the vortex. However, the creation of 
multiquantum vortices is impossible just by rotating the 
harmonic trapping potential at a high frequency, since 
the existence of many singly quantized vortices is energet- 
ically more favorable than a single multiquantum vortex 
0- Indeed, vortex lattices composed of single-quantum 
vortices were observed in such experiments jj]. 

Verifying the proposal of topological phase engineering 
by Nakahara et al. Q, Leanhardt et al. have recently 
succeeded in creating vortices simply by reversing the 
bias magnetic field used to trap the condensate. The 
vortices created display winding numbers with n = 2 (47r 
phase winding) or with n = 4 (8n phase winding), de- 
pending on the hyperfine spin states used for 23 Na con- 
densates. It was confirmed that the axial angular mo- 
mentum per particle is 2h (Ah) for the doubly (quarti- 
cally) quantized vortex. 

In this experiment, after creating a vortex it is held for 
~ 20 ms, watching the vortex core to split into single- 
quantum vortices. However, no splitting was observed 
in this time span. Since a doubly quantized vortex is 
expected to decay spontaneously into two singly quan- 
tized vortices owing to energetics, this observation seems 
to be puzzling. This motivates our investigation of the 
detailed dynamics of the decay process of multiply quan- 
tized vortices in view of the present experimental situ- 
ation. In dilute BECs, these were the first experimen- 
tally realized double-quantum vortices which are known 
to exist in other superfluids, such as superfluid 3 He-A 0. 
Therefore, we have a unique opportunity to examine the 
physics of multiply quantized vortices. 

There are several theoretical investigations on the in- 



stability of the doubly quantized vortex: The instability 
due to the bound state in the vortex core was pointed 
out by Rokhsar Q , who claimed that the decay of vortex 
states requires the presence of thermal atoms. On the 
other hand, Pu et al. @ found that even in the absence 
of thermal atoms, i.e. in the non-dissipative system, ap- 
pearance of the modes with complex eigenvalues leads to 
the spontaneous decay of multiply quantized vortices. 

Here we first study the Bogoliubov eigenvalue problem 
for a doubly quantized vortex in a cylindrically symmet- 
ric system in which the existence of modes with com- 
plex eigenfrequency, i.e., spectral instability, depends on 
the interaction strength. Then we investigate full three- 
dimensional cigar-shaped systems similar to the MIT ex- 
periments || . It is found that the splitting of the vortex 
core is initiated in different positions along the vortex 
axis, resulting in entangling vortices. 

Let us consider the spectral stability problem of the 
doubly quantized vortex by finding the spectrum of the 
quasiparticle excitations of the vortex state described by 
the order parameter ^(r). We apply the Bogoliubov 
eigenvalue equations [t| 

(C-fiLo a -/i)u Q (r)+ 5 [$(r)]V(r) = 0, 
#*(r)fu a (r) + (£ + H-,i)t. a (r) = 0, (1) 

where u> a are the eigenfrequencies related to the normal- 
mode functions u a (r) and v a (r). The chemical potential 
for the system of N bosons each having mass m is denoted 
as [i and the interaction strength is described by g = 
4Trh 2 a/m, where the s-wave scattering length is denoted 
as a. The operator £ = -h 2 V 2 /2m + V tI (r) + 2g|$(r)| 2 
consists of the kinetic energy, the trapping potential and 
the interaction term. In this paper, the trapping poten- 
tial is of the form Vtr(r) = m(uj 2 x 2 + uj 2 y 2 + lu 2 z 2 )/2, 
where u>i is the frequency of the trap in direction i. The 
doubly quantized vortex state created in the experiments 
is obtained by finding the minimum-energy solution 
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of the stationary version of the time-dependent Gross- 
Pitaevskii (GP) equation 
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with the constraint f|$(r, t)| <fr = N and ansatz 
5>(r,i) = </>(r, z)e z2e [lj], where (r,9,z) denote cylin- 
drical coordinates. The ansatz for the order parame- 
ter leads to constraints for the normal-mode functions 
u a {r) = w Q (r,z)e^+ 2 ) e , v a (r) = v a (r, z)e i( - q -^ e . Since 
the Bogoliubov eigenvalue equation is non-Hermitian, the 
existence of complex eigenvalues is not ruled out. If such 
an eigenvalue exists, the state is said to be spectrally 
unstable and small perturbations of the order parameter 
may grow exponentially in time. 

The Bogoliubov eigenvalue spectrum for double- 
quantum vortex states is computed for cylindrically sym- 
metric condensates with no phase modulation along the 
z-direction using a harmonic trapping potential u> x — 
Lu y = uj±_ and uj z = 0. In this case, there exist complex- 
frequency eigenmodes with q = ±2 for certain values 
of the interaction strength. The imaginary part of the 
eigenfrequency is plotted in Fig.QJa) as a function of the 
interaction strength an z = a J \4>(r)\ 2 dx dy . 

The emerging complex eigenvalues imply that by 
adding an infinitesimal perturbation, a pair of modes 
with complex conjugate eigenvalues is created sponta- 
neously and grow exponentially in time. In numerical 
calculations trap anisotropy A is used as the perturba- 
tion. One of these modes has q = — 2 and is confined 
into the vortex core. This mode describes the instabil- 
ity of the double-quantum vortex and leads to its split- 
ting into a pair of single-quantum vortices Q. The re- 
maining mode with q = +2 is created in order to ensure 
energy conservation. When the strength of interactions 
is weak (an z < 3), this mode corresponds to the co- 
rotating quadrupole mode. As the interaction strength 
an z increases, however, this mode changes from a node- 
less mode to the nodal mode along the r-axis and the 
real part of its frequency grows. 

To verify the role of the complex modes, using a 
slightly anisotropic trap with oj v = uj x / A — 2ir x 250 Hz 
(A = 1.01), we first solved the time-dependent GP equa- 
tion Pjl. which contains no term which can describe dis- 
sipative effects. The time evolution of the density pro- 
files with an z = 1.5 is shown in the upper panels of 
Fig. ^b). At t ~ 20 ms, the vortex splits into a pair 
of single-quantum vortices and the condensate simulta- 
neously changes into quadrupolar shape. Due to the co- 
rotating quadrupole excitation with q = +2, the conden- 
sate rotates with a constant frequency and "ghost" vor- 
tices ^3 nucleate near the condensate surface, which is 
clearly seen through the density modulation around the 
surface at t = 36 ms. Owing to the absence of dissipa- 
tion, the pair of single-quantum vortices precesses around 
the trap center and cannot escape from the trap. 
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FIG. 1: (a) Imaginary part of the complex excitation fre- 
quency in cylindrically symmetric ground states as a func- 
tion of an z . (b) Time-evolution of the density profiles with 
an z = 1.5 (upper panels) and 13.75 (lower panels) in the ab- 
sence of dissipation, (c) In the dissipative system (7 = 0.005), 
the time-development of angular momenta for several parti- 
cle numbers, corresponding to (1) — (6) in (a). Insets in (c) 
represent contour plots of particle density with an z — 13.75. 



For the interaction strength corresponding to the ex- 
periments (an z — 13.75), however, it is seen that within 
the time interval where the condensates are held ex- 
perimentally the splitting does not occur clearly in this 
cylindrically symmetric case shown in the lower panels 
of Fig. ^b). It should be noted that the characteristic 
time for the splitting is proportional to the inverse of the 
imaginary part of the eigenfrequency, but the mode ex- 
cited with q = — 2 has not yet split the vortex in Fig.QJb). 
The above calculations have also been carried out for the 
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anisotropics A = 1.02 and 1.03, but no difference was ob- 
served in the splitting of the vortex. The detailed depen- 
dence of the splitting on trap anisotropy is left for future 
research. 

In order to investigate the dissipative dynamics, we 
introduce a phenomenological damping term [lOl 1 1 1| to 
the GP equation through the substitution t — ► (1 — i"f)t, 
7 = 0.005. The temporal development of the angular mo- 
mentum is plotted in Fig-Q^c) f° r several values of the in- 
teraction strength and dynamics of the double-quantum 
vortex with an z — 13.75 is shown in the insets. It is 
found that when modes with a complex eigenvalue exist, 
the angular momentum decreases with time, which is re- 
lated to the decay into a pair of single-quantum vortices 
which spiral out to the condensate surface. In contrast, 
when there occur no complex eigenvalues, the double- 
quantum vortex becomes stable for a longer time. 

The experimentally realized multiquantum vortices 
were created in a cigar-shaped condensate confined before 
the vortex-generation process in a harmonic trap with the 
frequencies uj± = 2ir x 210 Hz and uj z = 2ir x 6 Hz |5j. 
After the creation, the radial trapping frequency was 
uj±_ = 2tt x 250 Hz and the condensate with N = 10 6 
atoms was unconfined along the z-direction. The uncon- 
finement is a technical difficulty which can be removed by 
reversing the magnetic field trapping the atoms along the 
z-direction when the bias field vanishes or by adding an 
optical potential along the z-direction so that the trap- 
ping frequency stays constant. Since it is possible to re- 
alize non-expanding condensate with a double-quantum 
vortex and the Bogoliubov equations can not be solved 
for an unconfined case, we restrict our computation to the 
confined case with ui± = 2tt x 250 Hz and lo z — 2tt x 6 Hz. 
Even though the condensate expands in the experiments 
, the qualitative behavior of the splitting of the double- 
quantum vortex is revealed by our results. 

After finding the minimum-energy solution for the con- 
densate with a double-quantum vortex along the z-axis 
|l4|. we plotted the strength of the interactions an z (z) 
ii J \tfi(z,x,y)\ 2 dxdy as a function of the z-coordinate in 
Fig. GJa) , from which it is seen that the linear densities 
n z in the regions \z\ < 51 /im and \z\ > 142 /an cor- 
respond to the two leftmost complex eigenfrequency re- 
gions in Fig.^a). Therefore, it is expected that the split- 
ting of the vortex begins in those regions and that the 
Bogoliubov spectrum features many complex-frequency 
modes, since the cigar-shaped condensate may be visu- 
alized to consist of several non-interacting condensates 
in icy-planes if one neglects the kinetic energy along the 
z-direction. 

The above argument is verified by first solving the Bo- 
goliubov equations for the cigar-shaped rotational BEC, 
which show that the vortex state features many complex- 
frequency eigenmodes with q = —2, plotted in Fig. 03a). 
The frequencies possess noticeable finite imaginary parts 
when their real part is in either one of two distinct do- 
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FIG. 2: (a) Linear density n z = J \cj>(x,y, z)\ 2 dx dy at 
t — 0. The vertical axis corresponds to the horizontal axis 
of Fig. ^a) . (b) Isosurfaces of the particle density and the 
average particle densities in the z-direction n±(x,y) demon- 
strating the splitting of the doubly quantized vortex at t — 
18.1 and 41.4 ms with anisotropy A = 1.01. 
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FIG. 3: (a) Eigenfrequencies for the cigar-shaped condensate 
having q = —2 plotted in the complex plane. Contour plots 
of density distribution of the modes marked with a rectangle 
and circle are plotted in subfigures (b) and (c) respectively. 



mains which are related to the two regions in the z- 
coordinate where the splitting of the vortex is expected to 
begin. The mode marked with a rectangle in Fig.^a), is 
localized in the region \z\ < 60 /im (see Fig. EHb)), while 
the other complex modes with lm(uj a ) > 0.02 occupy the 
region |z| > 100 ^m (see Fig.^c)). 

We solve the non-dissipative temporal evolution of the 
doubly quantized vortex state using the GP equation 
with the trap anisotropy A = 1.01. Isosurface plots 
of the particle density after 18.1 and 41.4 ms of time 
propagation are shown in Fig.[2Ib), from which it is ob- 
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served that the doubly quantized vortex had split into 
two single-quantum vortices first at the ends and then in 
the middle of the condensate, consistent with the above 
discussion. Also the particle density integrated in the 
^-direction n±(x, y) — J \4>{x, y, z)\ 2 dz is displayed to il- 
lustrate the density distribution observed in the experi- 
ments 0. If a dissipation term 7 = 0.005 is used, the 
splitting is slightly faster, but the shape of the condensate 
is not deformed in time as much as in the non-dissipative 
system, since the surface waves seen at the ends of the 
condensate are damped. 

Simulations performed for the cylindrically symmet- 
ric system show that an off-axis vortex tends to precess 
in the positive direction and the frequency of precession 
is an increasing function of the precession radius [l2j . 
However, numerical simulations of the GP equation show 
that the repulsion between the two vortices in a vortex 
pair renders the frequency of precession to be a decreas- 
ing function of the precession radius. Analogously to the 
cylindrically symmetric system, two vortices tend to twist 
around each other in a cigar-shaped BEC such that the 
parts where the distance between the vortices is large lag 
in precession compared with the parts where the distance 
is small, resulting in an entangled vortex pair [lfij . 

Following the above argument, the twisting of the vor- 
tices at t = 41.4 ms is negative in the region 75 //m < 
z < 150 /im, i.e., the azimuthal angle of the position of 
the vortices is a decreasing function of the z-coordinatc. 
The azimuthal angle decreases approximately by 107r as 
the z-coordinate changes from z — 78 /im to 130 /zm. 

Even though the vortex had clearly split in the mid- 
dle and at the ends of the condensate, only a single 
vortex core is seen in the z-integrated particle density. 
The reasons are that the twisting of the two singly quan- 
tized vortices causes the splitting to average out in the 
z-integration and that the vortex remains doubly quan- 
tized in the range 60 /im < z < 120 fim which covers a 
significant fraction of the condensed particles. This sug- 
gests an explanation of why no splitting was observed in 
the experiments |5j. 

In conclusion, the stability and the splitting of the 
double-quantum vortex created in Ref. Q has been an- 
alyzed via the Bogoliubov excitation spectrum and the 
time-dependent GP equation. The vortex is found to be 
spectrally unstable and the splitting of the vortex is in 
close analogy with the calculations for cylindrically sym- 
metric BECs. The emerging two single-quantum vortex 
lines tend to twist around each other, which makes it 
difficult to study the splitting of the doubly quantized 
vortex via particle density graphs integrated in the di- 
rection of the vortex line. Therefore, we suggest particle 
density profiles integrated in a direction perpendicular to 
the vortex line to be extracted [l^j ] in future experiments. 
We also encourage an experiment with different total par- 
ticle numbers to investigate the various combinations of 
instability regions shown in Fig. ^a). Our simulations 



show that the splitting of the doubly quantized vortex 
begins only from the ends of the condensate for peak in- 
teraction strength an z < 3. A confining potential should 
also be used as suggested to increase the lifetime of the 
condensate. 
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